/*******************************************************************************

        SSC0143 - Programacao Concorrente
        Prof. Julio Cesar Estrela
       
        T1 - Calculo de PI com 10000000 de casas decimais
       
        Grupo:
       
        Daniel de Oliveira Pereira              6427171
        Guilherme Takasawa Yagui                6426698
        Vinicius Marques Stocco                 6427233
       
        Março de 2013
       
********************************************************************************        
                 Algoritmo de Gauss-Lagendre - Versao Sequencial - tempo = 80.130s       
*******************************************************************************/        
        
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<gmp.h>

// numero de casas decimais corretas de PI
#define PREC_DGT        10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas
#define PREC_BIT        (int) ceil((log(10) / log(2)) * PREC_DGT)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado
#define ITR             (int) ceil(log(PREC_DGT) / log(2))


int main(void)
{	
	mpf_t a, b, t, p, pi, aux;
	int i = 0;
	mpf_set_default_prec(PREC_BIT);

        // marca  o clock inicial da execucao do algoritmo
	clock_t tempo = clock();

	mpf_init(a);
	mpf_init(b);
	mpf_init(t);
	mpf_init(p);
	mpf_init(aux);
	mpf_init(pi);

	//a(0) = 1
	mpf_set_d(a, 1.0);

	//b(0) = 1/sqrt(2)
	mpf_set_d(b, 1.0);
	mpf_set_d(aux, 2.0);
	mpf_sqrt(aux, aux);
	mpf_div(b, b, aux);

	//t(0) = 1/4
	mpf_set_d(t, 1.0);
	mpf_div_ui(t, t, 4);

	//p(0) = 1
	mpf_set_d(p, 1.0);
	
	while(i < ITR)
	{
		//pi = ((a(n)+b(n))^2)/(4*t(n))
		mpf_add(pi, a, b);
		mpf_pow_ui(pi, pi, 2);
		mpf_div(pi, pi, t);
		mpf_div_ui(pi, pi, 4);

		//gmp_printf("Iteracao %d, pi = %.*Ff \n", i, PREC_DGT, pi);

		mpf_set(aux, a);

		//a(n+1) = (a(n)+b(n))/2
		mpf_add(a, a, b);
		mpf_div_ui(a, a, 2);

		//b(n+1) = sqrt(a(n)*b(n))
		mpf_mul(b, aux, b);
		mpf_sqrt(b, b);

		//t(n+1) = t(n)-p(n)*(a(n)-a(n+1))^2
		mpf_sub(aux ,aux ,a);
		mpf_pow_ui(aux ,aux ,2);
		mpf_mul(aux ,p ,aux);
		mpf_sub(t, t, aux);

		//p(n+1) = p(n)*2
		mpf_mul_ui(p, p, 2);

		i++;
	}

	mpf_clear(a);
	mpf_clear(b);
	mpf_clear(t);
	mpf_clear(p);
	mpf_clear(aux);
	mpf_clear(pi);

	// calcula o numero clocks necessario para executar o algoritmo
        tempo = clock() - tempo;

	printf("\nTempo de execucao : %.3f s \n\n", (double) tempo / CLOCKS_PER_SEC);
	
	return 0;
}
